PROGRAM c3dFoam
use Acoustics, only: AcousticSurf, nNoiseSignals,NoiseSignal
use Data_cfoam
use DataTime
use Data_Mesh
use Data_Transducer !, only: defCellsTransducer
use DataCells
use DataFaces
use DataPerfectGas, only: Cv, nu, BulkVisc, BulkViscTotal, DynamicViscosity, TotalBulkViscosity
!use signalC, only: flag_dump
use TimeMeasure, only: GetHWTimerId,GetHWTick,StartHWTimer,StopHWTimer
use DataViscosity, only: CalcViscosityOn
use DataSource
use version
use InletBoundaryConditions
implicit none
    include 'mpif.h'

    character(256) params(2),starttime,finishtime
    character(1024) :: buf

    real *8 p00,tout,tin,localRe,rr,rr1,rr2,xt(3,8) ! ,defFactor
    real(8) xLocation(3) ! coordinates of observer
    real(8) dUdt(3),drdt,dpdt

    integer jpropeller,i1,ndeltap,j1 !,lperiodic(4) !,ntimemax ! for rotating coorsystem
    real *8 sumFX1,sumFX2 !, &! for postprocessing
    real *8 sumMX1,sumMX2 ! for postprocessing
    real *8 force(3) ! for rotating coorsystem
    real *8 rbuf(20000)  ! for postprocessing
    real *8 array(4)! for postprocessing
    real *8     cZ(3) ! coord of cell centers
    real *8 Sval(3)
    real *8 vcoriolis(3) ! for rotating coorsystem

    real(8) :: r(3), c, time, Q,dQdT,P0,U0(3), Pi,ek,a, Rho0, timeMem, timeMem1

!  c - center
!  f - forward (next cell)
!  s - surface
!  sf - middle time step

	real *8 xg(3) ! local node's coord ? maybe face's centers
	real *8 grad(3,3) ! local tensor diff velocities Dv_i/Dx_j
	real(8) :: DynVisc ! dynamic viscosity
	real(8) :: BulkViscousStress ! bulk viscous stress
	real(8) :: DivVStress ! divergence (Velocity * Viscous_Stress)
!	real *8 nu ! local coef viscosity
!	real(8) :: nu1,nu2 ! auxilary viscosities for smoothed simulations
	real *8 stress(3,3) ! local Reynold's stess
!	real(8) :: divU ! divergence part in viscid stress tensor
	real *8 uf(3), ucf(3)  ! local velocities' basis
	real *8 ucsf(3) ! old value of ucf ! s = staroe!
	
	real(8) :: memR
    real *8,pointer::cellOut(:,:)
	character *4 char48
	
	integer i,j,k,l,m !,IERR,rank
	real *8 p1(2,100)
	
    real *8 Dbmin,Dfmin ! corrector to current and next face
    real *8 delta,distance ! geometrical factors
    real *8 ds1,ds2,dsmall !,dt ! for determining time step
    real *8 dt2p,dtmin,dtminx,dtn ! dt2 for determining time step ???
    real *8 dtp ! dtOut for output
!    real *8 epsilon ! numerical dumping (viscosity)
    real *8 factor  ! isn't used
    real *8 ttimeOut   ! for output
    real *8 uin ! for inflow
    real *8 vmax ! max velocity
    real *8 vol ! volume
    real *8 vtmp ! mas  flux, Rho*U
    real *8 xg1  ! node's coord for postprocessing
    real *8 xg5  ! node's coord for postprocessing
    real *8 omega  ! rotate velocity for rotating system
    ! MPI
    integer ier,irank,irank1,is,isdispls,jj,irank0,jstart,lsmall,& ! ie,
        m1,ndelta,netr0,nodeInterfacet,& ! ,nCellsOut,netp0,nodeInterface
        nout,np1,iconvname,n,nOut1
    integer :: iProc,statusp(mpi_status_size) 
!    integer :: nUnusedSections ! number of sections with more than 4 points, for acoustics surface
    integer :: nCellsOutSum !, nUnusedSectionsSum ! for parallel
	real *8 mu ! viscosity, not used
    real *8 khi ! cinductivity, not used
    real *8 accel ! for rotating system
    real *8 rwall ! for postprocessing
    real *8 xlc(3,10) ! local vars for postprocessing
	real *8 xc(3),scale !,sound ! local vars for postprocessing

    ! output
	integer  nstop  ! maximum of numbers of time steps
	integer  ivertex   ! ?

    real(8) :: alfa, NoSlipCoef
    
    integer :: iCell,iFace ! indeces for cycle on Cells and Faces

!	data lperiodic/1,4,3,2/
	
    character *3 :: DirString
    integer :: Counter, kk, nPrintOut

    real(8) :: auxReadReal, aa

    real(8) :: UniformNormP,UniformNormU, L1NormP,L1NormU, L2NormP,L2NormU, L2AveNormP,L2AveNormU
    real(8) :: deltaP, deltaU, Volume, Umax, Pmax

    integer :: nStepCellsMin,nStepCell, nStepCellsMin1, nCellsSum ! ,nStepCellsSum
    real(8) :: nStepCells1, nStepCellsSum, maxDbCell,maxDbCell1, minDbCell,minDbCell1, nStepHomogeneous
    integer :: iCellFace
    
    integer :: bmeshs, bmeshf

! Timers
    integer :: Timer_ReadData
    integer :: Timer_TransducersInit
    integer :: Timer_ReadMesh
    integer :: Timer_MeshOut
    integer :: Timer_SetInitCond
    
    integer :: Timer_OutRes
    integer :: Timer_Phase1
    integer :: Timer_Phase2
    integer :: Timer_Phase3
    integer :: Timer_OutputTransducers
    integer :: Timer_Acoustics
    
    real :: Secs
	integer :: Ncalls
	integer(kind=8) :: nfreq
	character(16) :: TimerName

    integer :: OpenFile
    
    real(8),allocatable :: ViscousFluxes_view(:,:,:), faceGradU_view(:,:,:) ! sumViscousCellFluxes(:,:),

    integer(8) time1, time2, timebase

!DEC$ IF DEFINED (__INTEL_COMPILER)
	call cpu_time(CPU_Start_Time)
!DEC$ ELSE
    CPU_Start_Time = mpi_wtime()
!DEC$ ENDIF

	! MPI initialization
    call MPI_INIT(IERR)
    if(ierr.ne.MPI_SUCCESS)stop 'MPI initilization failed'
	call MPI_COMM_SIZE ( MPI_COMM_WORLD, nump, ierr)
    if(ierr.ne.MPI_SUCCESS)stop 'MPI initilization failed'
	CALL MPI_COMM_RANK ( MPI_COMM_WORLD, rank, ierr)
    if(ierr.ne.MPI_SUCCESS)stop 'MPI initilization failed'

    call GetArg(1,buf)
    if(buf(1:1).ne.' ')then
        fname = trim(buf)
    end if
        
    if(rank.eq.0)then
    
        call WriteProgramVersion()
        
        call Date_and_Time(params(1),params(2))
        starttime=params(1)(7:8)//'.'//params(1)(5:6)//'.'//params(1)(1:4)   &
        //' at '//params(2)(1:2)//':'//params(2)(3:4)//':'//params(2)(5:6)

        write(*,*) 'Simulation start '//trim(starttime)
        write(*,*) 'Input file    '//trim(fname)
        write(*,*)

        call date_and_time(TIME=Eta)
    end if

    call GetHWTimerId('ReadData',Timer_ReadData)
    call GetHWTimerId('TransducersInit',Timer_TransducersInit)
    call GetHWTimerId('ReadMesh',Timer_ReadMesh)
	call GetHWTimerId('MeshOut',Timer_MeshOut)
	call GetHWTimerId('SetInitCond',Timer_SetInitCond)

	call GetHWTimerId('OutResults',Timer_OutRes)
    call GetHWTimerId('Phase1',Timer_Phase1)
    call GetHWTimerId('Phase2',Timer_Phase2)
    call GetHWTimerId('Phase3',Timer_Phase3)
    call GetHWTimerId('OutputTransducers',Timer_OutputTransducers)
    call GetHWTimerId('Acoustics',Timer_Acoustics)

    ! for acoustics postprocessor
    Pobs=0d0
    xLocation=0d0
    xLocation(2)=30*0.05d0 !28.
        
    utmp(:,:) = 0.0
    
	! Preprocessing
	nout=0
	btype=0

	call StartHWTimer(Timer_ReadData)

	call ReadDataFile(gam,mu,khi,CFL,UFON,VFON,WFON,ROFON,PFON,TFON, &
                         nstop,nprint,btype,utmp,ivertex,xc,irestart,scale,xlc,accel &
                         ,sound,epsilon,rwall,omega,dtOut, &
                         ttime,dt,TimeStepMax,TimeStepMin)

    call CalcNozzleBoundCondInit()

	call StopHWTimer(Timer_ReadData)

	dt=TimeStepMin ! initial time step
	mu = mu / rofon
	BulkVisc = BulkVisc/ rofon
	nu=mu
	if(nu.gt.ZeroEpsilon)then
	    BulkViscTotal = (BulkVisc - nu * 2d0 / 3d0) * 0.5 / nu ! multiplicator 0.5 / nu because bulk viscosity added to Grad_i U_j
	                                                       ! where gradU(i,j,k)=nu*(grad(i,j)+grad(j,i))
    else
        BulkViscTotal = 0d0
    end if

    if(AcousticsOnOff)then
	    if (rank.eq.0) then ! Georgiy's new acoustics postprocessor
                allocate (Pobsw(300000*nump))
        else
                allocate (Pobsw(1))
        end if
    end if

	allocate(sendCounts(nump+1))
	allocate(sDispls(nump+1))    
	allocate(recvCounts(nump+1)) 
	allocate(rDispls(nump+1))    
	allocate(sendCountsv(nump+1))
	allocate(sDisplsv(nump+1))   
	allocate(recvCountsv(nump+1))
	allocate(rDisplsv(nump+1)) 
	allocate(sendCountsVPoints(nump+1)) ! for sending VPoints
	allocate(recvCountsVPoints(nump+1)) 
	allocate(sDisplsVPoints(nump+1)) 
	allocate(rDisplsVPoints(nump+1))
	allocate(sendCountsVPoints9(nump+1)) ! for sending VPoints
	allocate(recvCountsVPoints9(nump+1)) 
	allocate(sDisplsVPoints9(nump+1)) 
	allocate(rDisplsVPoints9(nump+1))
	allocate (BoundaryProcessorNumber(nump+1))  
    sendCounts=0
    sDispls=0
    recvCounts=0
    rDispls=0
    sendCountsv=0
    sDisplsv=0
    recvCountsv=0
    rDisplsv=0

	write(char48,'(i4)') rank
	is=4
	if (rank.gt.9) is=3	
	if (rank.gt.99) is=2	
	if (rank.gt.999) is=1	
	i=5-is
!	processor(9+1:9+i)=char48(is:4)
	processor='processor'//char48(is:4)//'/'
	ie=9+i+1
	if (nump.eq.1) then
	    processor='./'
	    ie=2
	end if
	! end MPI initialization
	! Read Boundary Mesh / Boundary Mesh has FOAM's format

	call StartHWTimer(Timer_ReadMesh)
	call ReadMeshOpenFOAM()
	call StopHWTimer(Timer_ReadMesh)

    call MPI_BARRIER(mpi_comm_world,ierr)
!    goto 20000

    if(CalcViscosityOn.and.AsynchronousTime)then
        call CalcPointsFacesCellsNeigForVisc()
    end if

    if(NewPeriodicNumeration)then
        call CalcNewPeriodicNumeration()
    end if
    
    if(SplitPeriodicBoundaries)then
        call PR_SplitPeriodicBoundaries()
    end if

	call StartHWTimer(Timer_TransducersInit)
    call defCellsTransducer() ! definition of cells for output data in Transducers
	call StopHWTimer(Timer_TransducersInit)

    outResults = 32
    outRestarts = outResults + 1

	call StartHWTimer(Timer_SetInitCond)
    call SetInitConditions()    ! end set initial conditions
    call StopHWTimer(Timer_SetInitCond)

!    call CalcJacobians()

	! initialization for MPI
	call StartHWTimer(Timer_MeshOut)
	open(1,file=processor(1:ie)//'outmeshhsp',form='unformatted')
!	write(*,*) rank,ncellsout
	write(1) nPoints,nCellsOut,((meshPoints(jj,i),jj=1,3),i=1,nPoints) &
                    ,((meshCellPoints(jj,meshCellOut(i)),jj=1,8),i=1,nCellsOut)
    close(1)
    call StopHWTimer(Timer_MeshOut)

	call StartHWTimer(Timer_OutRes)
	! Output files for results and restarts
    if(OutputFromOneProc)then
!        open(32,file='out',form='unformatted')
        if(rank.eq.0)then
            open(outResults,file='out',form='unformatted',IOSTAT=OpenFile)
!            open(32,file='out',form='unformatted',IOSTAT=OpenFile)
            if (OpenFile.ne.0) then
	            write(*,*) '@@ Fatal Error: File open error file out'
	            STOP '@@ Fatal Error!'
            end if
 	        open(outRestarts,file='outRes',form='unformatted')
 	    end if
 	    call OutputFromAllProcsInit()
    else
! 	    open(32,file=processor(1:ie)//'out',form='unformatted')
! 	    open(outRestarts,file=processor(1:ie)//'outRes',form='unformatted')
 	end if
	call StopHWTimer(Timer_OutRes)
	
    outFW = outRestarts + 1 ! nump +
    
 	open(outFW,file=processor(1:ie)//'outhsp',form='unformatted')
 	allocate (cellOut(5,nCellsOut))
 	cellOut=0.0d0
 	ttimeOut=dtOut
!	dt = 0.1d-7
!    dt = min(TimeStepMax,max(dt,TimeStepMin))
!	dt2=dt/2.0d0
    ndelta=400

    ! output MeshCellCenters and MeshFaceCenters
    if(JustMeshCellFaceCentersOut)then
        call WriteMeshCellFaceCenters()
        call MPI_Barrier(MPI_COMM_WORLD,ierr)
        Stop 'JustMeshCellFaceCentersOut'
    end if

    if(CalcBoundaryFaceCubesOnly)then
        call CalcBoundaryFaceCubes()
    end if

    if(OutputFaceCubesOnly)then
        call OutputFaceCubes()
    end if

	call StartHWTimer(Timer_OutRes)
    if(OutputFromOneProc)then
        call OutputFromAllProcs()
    end if
    call StopHWTimer(Timer_OutRes)

	call StartHWTimer(Timer_TransducersInit)
    call OutputTransducerDataInit()
    call StopHWTimer(Timer_TransducersInit)

	epsilonSmoothed = epsilon ! 1. !
	propor = 1.0d0

! Simulation start

    allocate(sumViscousCellFluxes(3,nCells))

    if(SourcePresent)then
        call CalcSourceLocation(SourceLocation, iCellSource)
    end if

! _________________________________________________________
! _____ complexity calculation ____________________________
! _________________________________________________________
    c = 340d0
    nStepCellsMin = 0
    timeMem = minVal(DbCells(1:nCells))

    if (nump.gt.1) then
        nCellsSum = 0
        call MPI_ALLREDUCE(timeMem,timeMem1,1,MPI_DOUBLE_PRECISION,MPI_MIN,MPI_COMM_WORLD,ierr)
        call MPI_Barrier(MPI_COMM_WORLD,ierr)
        timeMem = timeMem1

        call MPI_ALLREDUCE(nCells,nCellsSum,1,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,ierr)
        call MPI_Barrier(MPI_COMM_WORLD,ierr)

        maxDbCell = maxVal(Db(1:nCells))
        call MPI_ALLREDUCE(maxDbCell,maxDbCell1,1,MPI_DOUBLE_PRECISION,MPI_MIN,MPI_COMM_WORLD,ierr)
        call MPI_Barrier(MPI_COMM_WORLD,ierr)
        maxDbCell = maxDbCell1


        minDbCell = minVal(Db(1:nCells))
        call MPI_ALLREDUCE(minDbCell,minDbCell1,1,MPI_DOUBLE_PRECISION,MPI_MIN,MPI_COMM_WORLD,ierr)
        call MPI_Barrier(MPI_COMM_WORLD,ierr)
        minDbCell = minDbCell1
    else
        maxDbCell = maxVal(Db(1:nCells))
        minDbCell = minVal(Db(1:nCells))
        nCellsSum = nCells
    end if

    timeMem = 10d0 * maxDbCell ! maxVal(DbCells)
    nStepHomogeneous = timeMem / minDbCell ! minVal(DbCells) ! + 10

    nStepCells = 0
    nStepCells1 = 0d0
    do iCell = 1,nCells
        aa = int((timeMem+1d-12) / DbCells(iCell))
        if(timeMem.gt.aa*DbCells(iCell)) aa = aa + 1
        nStepCells1 = nStepCells1 + aa ! int((timeMem+1d-12) / DbCells(iCell)) !  + 1 ! / cfl + 1
    end do

    if (nump.gt.1) then
        nStepCellsSum = 0
        call MPI_ALLREDUCE(nStepCells1,nStepCellsSum,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
        nStepCells1 = nStepCellsSum
    end if

    if(rank.eq.0)then
        nStepCells = int((nStepCells1+1d-12)/nCellsSum)
        write(*,*)
        write(*,*) 'Theoretical maximum speed-up'
        write(*,'(a22,i15)') '    nStepsAsynchrony = ',nStepCells ! , ', rank=', rank
        write(*,'(a22,f15.0)') '    nStepsHomogenous = ',nStepHomogeneous ! , ', rank=', rank
        write(*,'(a22,f15.2)') '    maxAcceleration  = ',nStepHomogeneous/nStepCells ! , ', rank=', rank
        write(*,'(a18,f15.2,2f25.10)') '    maxDb/minDb = ', maxDbCell/minDbCell, maxDbCell,minDbCell ! , ', rank=', rank
        write(*,*)
    end if

    nStepHomogeneous = 1d0/minDbCell !/ nCellsSum
    
    nStepCells = 0
    nStepCells1 = 0d0
    do iCell = 1,nCells
        nStepCells1 = nStepCells1 + 1d0/DbCells(iCell) !/ nCellsSum
    end do
    nStepCells1 = nStepCells1 / nCellsSum
    
    if(rank.eq.0)then
        nStepCells = nStepCells1
        write(*,*)
        write(*,*) '***************************************************************'
        write(*,*)
        write(*,*) 'Theoretical maximum speed-up'
        write(*,'(a22,i15)') '    nStepsAsynchrony = ',nStepCells
        write(*,'(a22,f15.0)') '    nStepsHomogenous = ',nStepHomogeneous
        write(*,'(a22,f15.2)') '    maxAcceleration  = ',nStepHomogeneous/nStepCells1
        write(*,'(a18,f15.2,2f25.10)') '    maxDb/minDb = ', maxDbCell/minDbCell, maxDbCell,minDbCell
        write(*,*)
    end if

    if(nump.eq.1)then
        write(*,*) 'Theoretical maximum speed-up=', sum(DbCells(1:nCells))/(minDbCell*nCells)
    end if
    
    !stop 'Theoretical speedup'
    
! _____ complexity calculation finished ___________________________

    if(rank.eq.0)then
        do i = 1,5
		    call GetHWTick(i,Secs,Ncalls,freq=nfreq,name=TimerName)
		    if(ncalls.eq.0) cycle
		    write(buf,'(3a,i8,2(a,1pe12.3))') &
			    'Timer ',TimerName,' called ',NCalls,' times, CPU time =',secs, &
			    '; T/Ncalls =',secs/NCalls
		    write(*,*) trim(buf)
        end do
    end if

    do i=1,nCells
        pnew(i)=p(i)
        Tnew(i)=T(i)
        Unew(1:3,i)=U(1:3,i)	
        Enew(i)=E(i)          ! 
        rhonew(i)=rho(i)      !!
    end do

    call MinMaxWrite()

    if(writeInitCondOn)then
        if(OutputFromOneProc)then
            call OutputFromAllProcs()
        else
           call WriteRes()
        end if
    end if

!DEC$ IF DEFINED (__INTEL_COMPILER)
	call cpu_time(CPU_Start_Steps_Time)
!DEC$ ELSE
    CPU_Start_Steps_Time = mpi_wtime()
!DEC$ ENDIF

    if(AsynchronousTime)then
        call OutputTransducerAsynData()
    end if
    
    if(PreprOnly)then
    
        if(rank.eq.0)then
            call Date_and_Time(params(1),params(2))
            finishtime=params(1)(7:8)//'.'//params(1)(5:6)//'.'//params(1)(1:4)   &
            //' at '//params(2)(1:2)//':'//params(2)(3:4)//':'//params(2)(5:6)

            write(*,*) 'Simulation finish '//trim(finishtime)

!DEC$ IF DEFINED (__INTEL_COMPILER)
	        call cpu_time(CPU_Finish_Time)
!DEC$ ELSE
            CPU_Finish_Time = mpi_wtime()
!DEC$ ENDIF
            cpu_dtime = CPU_Finish_Time-CPU_Start_Time
	        write(*,*) 'CPU time=', cpu_dtime, ' sec'
	        write(*,*)

            write(*,*) 'Pre-processing done sucessfully!'
            write(*,*)
        end if

	    call MPI_FINALIZE(IERR)

        Stop
    end if

    do nstep=nStepInitial,nstop

	    if(AsynchronousTime)then
!	        ttime = 0.05d0 ! ttime + 2500d0 * dtOut	        
!	        ttime = ttime + dtOut
!	        ttime = 1.d0
            nDtGlobal = int((dtOut+TimeEpsilon)/TimeStepMin)
	        call CalcLocalSteps() ! forget about Acoustics and Transducers
	        ttime = ttime + dtOut

            if (mod(nstep,nprint).ne.0) cycle

            if (rank.eq.0)then
                write(*,'(a21,e12.5)') 'Cubes output, ttime =',ttime
!DEC$ IF DEFINED (__INTEL_COMPILER)
	                call cpu_time(cpu_2)
!DEC$ ELSE
	                cpu_2 = mpi_wtime()
!DEC$ ENDIF
	                cpu_dtime = cpu_2-CPU_Start_Steps_Time
	                write(*,*)'CPU time=', cpu_dtime
            end if

!            cycle

! ************************************************************* !
            do iCell = 1,nCells
                Tnew(iCell) = pnew(iCell)/(rhonew(iCell)*cv*(Gam-1.0d0))
            end do
            do iFace = 1,nFaces
                psnew(iFace) = FaceVarsNew(5,1,iFace)
                Usnew(1:3,iFace) = FaceVarsNew(2:4,1,iFace)
                Rhosnew(iFace) = FaceVarsNew(1,1,iFace)
            end do
! ************************************************************* !

            call MinMaxWrite()

            if(OutputFromOneProc)then
                call OutputFromAllProcs()
            else
                call WriteRes()
            end if
	        
	        cycle
	        
	    else
	        dtmin = TimeStepMax / cfl

            ! StepCalculation
		    if(nstep.lt.nTimeStepAcceleration) then ! ntimemax
			    propor = (nstep-1.0)/(nTimeStepAcceleration-1.0) ! ntimemax
			    epsilonSmoothed = 1.0d0-propor+propor*epsilon ! utmp(6,100)
		    end if

	        do k=1,nFaces
                Ps(k)=Psnew(k)
                Ts(k)=Tsnew(k)
                Us(1:3,k)=Usnew(1:3,k)
                Es(k)=Esnew(k)       ! 
                Rhos(k)=Rhosnew(k)  !!
                
	        end do
    ! phase 1
    
	        call StartHWTimer(Timer_Phase1)
            include 'phase1_sak.H' ! conservative substep
            call StopHWTimer(Timer_Phase1)

	        call StartHWTimer(Timer_OutputTransducers)
            call OutputTransducerData()
            call StopHWTimer(Timer_OutputTransducers)

    ! phase2 ! flux spliting
	        call StartHWTimer(Timer_Phase2)
            call CalcFacesVariables()
            call StopHWTimer(Timer_Phase2)

!	        do k=1,nFaces
!                if(isnan(Rhos(k)).or.isnan(Ps(k)).or.isnan(Us(1,k)).or.isnan(Es(k)))then
!                    write(112,*) 'alarm: Ln 732', nstep,Rhos(k),Ps(k),Us(1,k),Es(k)
!                    stop 'Nans are here 2'
!                end if
!            end do

    ! phase 3 ! Corrector
	        do i=1,nCells
                p(i)=pnew(i)
                T(i)=Tnew(i)
                U(1:3,i)=Unew(1:3,i)	
                E(i)=Enew(i)          ! 
                rho(i)=rhonew(i)      !!
	        end do

!call MinMaxWrite()
!call WriteRes()

	        call StartHWTimer(Timer_Phase3)
            include 'phase1_sak.H'
	        call StopHWTimer(Timer_Phase3)

        ! Output

!            call MinMaxWrite()

            if(AcousticsOnOff) then
                call StartHWTimer(Timer_Acoustics)
                if(AcousticSurf)then
                    call CalcNoiseSignal()
                else
                    include 'outFW_Georgiy.H' !'outFW.H' ! postprocessing ! Grigoriy's new acoustic postprocessor
                end if
                call StopHWTimer(Timer_Acoustics)
            end if

	        call StartHWTimer(Timer_OutputTransducers)
            call OutputTransducerData()
	        call StopHWTimer(Timer_OutputTransducers)

	        do i=1,nCells
                p(i)=pnew(i)
                T(i)=Tnew(i)
                U(1:3,i)=Unew(1:3,i)	
                E(i)=Enew(i)          ! 
                rho(i)=rhonew(i)      !!
	        end do

	        if (mod(nstep,nprint).eq.0) then

                call MinMaxWrite()
            
                call StartHWTimer(Timer_OutRes)
                if(OutputFromOneProc)then
                    call OutputFromAllProcs()
                else
                    call WriteRes()
                end if
                call StopHWTimer(Timer_OutRes)

	            if (rank.eq.0)then
	                write(*,'(a18,3e12.5)') 'ttime,dt,dt2=',ttime,dt,dt2

!DEC$ IF DEFINED (__INTEL_COMPILER)
	                call cpu_time(cpu_2)
!DEC$ ELSE
	                cpu_2 = mpi_wtime()
!DEC$ ENDIF
	                cpu_dtime = cpu_2-CPU_Start_Steps_Time
	                write(*,*)'CPU time=', cpu_dtime
	            end if
	        end if

        end if

    	if (mod(nstep,100).eq.0.and.rank.eq.0) then
            call date_and_time(TIME=Eta)
	        write(*,'(2i8,2e12.5,a20)') nstep,rank,ttime,dt,eta
	    end if

	    if(AsynchronousTime.eqv..true.) cycle
	    
!	    if(AsynchronousTime.eqv..false.)then
            ! Time Step Calculation
            if (nump.gt.1) then
	            call MPI_Barrier(MPI_COMM_WORLD,ierr)
                call MPI_ALLREDUCE(dtmin,dtn,1,MPI_DOUBLE_PRECISION,MPI_MIN,MPI_COMM_WORLD,ierr)
	            dtmin=dtn
	        end if

            dt=dtmin*cfl ! max(dtmin*cfl, TimeStepMin) ! definition of time step
!            dt=min(dt,TimeStepMax)
            dt2=dt/2.0d0
!        end if
	    
!	    if(AsynchronousTime.eqv..false.)then
            if (dtmin.gt.1) then
	            write(*,*) 'dtmin=',dtmin, nstep
	            stop 'Fatal error in dtmin'
	        end if
!	    end if
    end do
! end of main cycle / Time Loop



! deallocate
    if(AcousticsOnOff)then
        deallocate(Pobsw)
    end if

    if(OutputFromOneProc)then
        if(rank.eq.0)then
            close(32)
            close(outRestarts)
        end if
    else
        close(32)
        close(outRestarts)
    end if
    close(outFW)

!	call inputoutput('0/p','1/p',p,nCells)
	
	if(rank.eq.0)then
	    write(*,*) ''
	    write(*,*) 'Calculations have been ended successfully!'
	    write(*,*) ''
	end if

	if(AsynchronousTime)then
	    deallocate(SubStepCellVars)
	    deallocate(CharsCellPresent)

        deallocate(CharsOld,CharsNew)
        deallocate(fluxFaces)
        deallocate(fluxFacesReserve)
        deallocate(fluxesFacePresent)
	end if

20000 continue
	call MPI_FINALIZE(IERR)
	
end PROGRAM c3dFoam
	
	function deff(w1,w1f,w1cf,w1csf,char3f,dt,de,epsilon)
	use Data_cfoam, only: epsMaximumPrinciple
	implicit none
	    real *8 deff
	    real *8 w1,w1f,w1cf,w1csf,char3f,dt,de,w1n
	    real *8 w1well,w1maxf,w1minf,w1cfold
	    real *8 epsilon
    !	include 'epsilon.h'
	    w1well=(w1csf-w1cf)/(0.5d0*dt)+char3f*(w1f-w1)/de
	    w1maxf=dmax1(w1,w1cf,w1f)+dt*w1well + epsMaximumPrinciple
	    w1minf=dmin1(w1,w1cf,w1f)+dt*w1well - epsMaximumPrinciple
    !	w1n=2*w1csf-w1f
	    w1n=(2*w1csf-w1f*(1.0d0-epsilon))/(1.0d0+epsilon)

	    if(w1n.gt.w1maxf) w1n=w1maxf
	    if(w1n.lt.w1minf) w1n=w1minf
	    deff=w1n
	    return
	end function deff
	    
    function defb(w1,w1b,w1cb,w1csb,char3b,dt,de,epsilon)
    use Data_cfoam, only: epsMaximumPrinciple
    implicit none
	    real *8 defb
	    real *8 w1,w1b,w1cb,w1csb,char3b,dt,de,w1n
	    real *8 w1well,w1maxb,w1minb,w1cbold
	    real *8 epsilon
    !	include 'epsilon.h'
	    w1well=(w1csb-w1cb)/(0.5*dt)+char3b*(w1-w1b)/de
	    w1maxb=dmax1(w1,w1cb,w1b)+dt*w1well + epsMaximumPrinciple
	    w1minb=dmin1(w1,w1cb,w1b)+dt*w1well - epsMaximumPrinciple
    !	w1n=2*w1csb-w1b
	    w1n=(2*w1csb-w1b*(1.0d0-epsilon))/(1.0d0+epsilon)
	    if(w1n.gt.w1maxb) w1n=w1maxb
	    if(w1n.lt.w1minb) w1n=w1minb
	    defb=w1n
	    return 
	end function defb
	
! ALCF: moved "coor" to a separate PR_coor.f90 file

	subroutine transform(a,b,c)
	implicit none
!	IMPLICIT REAL *8(a-h,o-z)
	    real(8),intent(in) :: a(3)
	    real(8),intent(out) :: b(3)
	    real(8),intent(in) :: c(3,3) 
	    integer j
	    do j=1,3
    !	b(j)=c(j,1)*a(1)+c(j,2)*a(2)+c(j,3)*a(3)
	        b(j)=c(1,j)*a(1)+c(2,j)*a(2)+c(3,j)*a(3)
	    end do
	    return
    end subroutine transform
    
    subroutine transformr(a,b,c)
    implicit none
    !	IMPLICIT REAL *8(a-h,o-z)
	    real *8 a(3)
	    real *8 b(3)
	    real *8 c(3,3)
	    integer j
	    do j=1,3
	        b(j)=c(j,1)*a(1)+c(j,2)*a(2)+c(j,3)*a(3)
    !    b(j)=c(1,j)*a(1)+c(2,j)*a(2)+c(3,j)*a(3)
	    end do
	    return
	end subroutine transformr

	subroutine inputoutput(name,name1,p,nump)
    implicit none
        character (*) name,name1
	    character *280 a
        real *8 p(1)
        integer nump
        integer i
        
	    open(11,file=name)
	    open(12,file=name1)
	    do i=1,19
	        read(11,'(a280)') a
	        write(12,'(a280)') a
	    end do
        write(12,*) 'internalField   nonuniform List<scalar>'
        write(12,*) nump
        write(12,*) '('
        do i=1,nump
            write(12,*) p(i)
        end do
        write(12,*) ');'
	    read(11,'(a280)') a
	    do i=1,32
	        read(11,'(a280)') a
	        write(12,'(a280)') a
	    end do
	    close(11)
	    close(12)
        return
    end subroutine inputoutput
    
	subroutine add_1(m,n1)	
    implicit none
	    integer m(8)
	    integer n1(4)
	    integer n
	    integer l(8)
        integer i,i1,j,i0,i01,j0,j00,j01,j02,ll
        
	    l=0
    !	write(*,*) 'n1=',n1
	    if (m(1).eq.0) then
	        m(1:4)=n1(1:4)
	        return
	    else
	        l=0
	        do i=1,4
	            do j=1,4
	                if (m(i).eq.n1(j)) then
	                    l(i)=1
	                    l(4+j)=1
	                end if
	            end do
	        end do
            ll=l(1)+l(2)+l(3)+l(4)
            if (ll.eq.0) return
            do i=1,4
                i1=i+1
                if (i1.gt.4) i1=1
                if (l(i)+l(i1).eq.2) then
                    i0=i
                    i01=i1
                end if
                if (l(i+4)+l(i1+4).eq.2) then
                    j0=i
                    j01=i1
                end if
            end do
            j00=j0-1
            if (j00.lt.1) j00=4
            j02=j01+1
            if (j02.gt.4) j02=1
	        if (m(i0).eq.n1(j0)) then
                m(i0+4)=n1(j00)
                m(i01+4)=n1(j02)
            else
                m(i0+4)=n1(j02)
                m(i01+4)=n1(j00)
            end if
        end if
	    return
	end subroutine add_1
	
	subroutine vikhr(xg,U,p,rho,etvikhr)
    implicit none
!	implicit none
	!IMPLICIT REAL *8(a-h,o-z)
        real *8 p
        real *8 rho
	    real *8 xg(3)
	    real *8 xc(3)
	    integer ivar
	    real *8 gam1
	    real *8 x
	    real *8 y
	    real *8 u(3)
	    real *8 rovikhr
	    real *8 pvikhr
	    real *8 etvikhr
	    real *8 uvikhr
	    real *8 vvikhr
	    real *8 eivikhr 
	    real *8 a
	    real*8 rofon1,pfon1,x0,y0,epsilon,alf,rc,rd,taut,sint,cost, &
    	    expo,dexpo,gamy,gamx,gam2
	    real *8 ufon,vfon,rofon,pfon,eifon,etfon
	    real *8 rr,v1,v2,gam,alfp,cv,dtmint
        real *8 deff,defb
            
	    gam=1.4d0
        gamy=((gam-1.0d0)/(2.0d0*gam))
	    gamx=2.0d0*dsqrt(gam)/(gam-1.0d0)
	    gam2=2.0d0*gam
        rofon1=1.0d0
        pfon1=1.0d0
	    x0=0.0d0
	    y0=0.0d0
	    x=xg(1)
	    y=xg(2)
	    epsilon=0.204d0
	    alf=0.3d0
	    rc=0.05d0
	    rd=dsqrt((x-x0)**2+(y-y0)**2)
	    taut=rd/rc
	    sint=(y-y0)/rd
	    cost=(x-x0)/rd
	    expo=dexp(alf*(1-taut**2))
	    dexpo=(gam-1.0d0)/4.0d0/alf/gam*(epsilon*expo)**2
	    rovikhr=rofon1*(1.0d0-dexpo)**(1.0d0/(gam-1.0d0))
	    uvikhr=epsilon*taut*expo*sint	
	    vvikhr=-epsilon*taut*expo*cost
	    pvikhr=pfon1*(1.0d0-dexpo)**(gam/(gam-1.0d0))
	    eivikhr=pvikhr/(gam-1.0d0)/rovikhr
	    etvikhr=eivikhr+0.5d0*(uvikhr**2+vvikhr**2)
	    u(1)=uvikhr
	    u(2)=vvikhr
	    u(3)=0.0d0
	    p=pvikhr
	    rho=rovikhr
	    return
	end subroutine vikhr
	
	subroutine sol_eq(a,rd,r0,alf,p)
	implicit none
	    real *8 a,rd,r0,alf,p
	    real *8 r,dr
	    integer n,i
	    
	    dr=1e-4
	    n=rd/dr+1
	    p=0.0d0
	    r=dr/2.0d0
	    do i=1,n
	        p=p+(a*r*dexp(-2*alf*r**2/r0**2))*dr
	        r=r+dr
	    end do
	    return
	end	subroutine sol_eq

	subroutine interp(p,n,z,rtek)
    implicit none
	    !IMPLICIT REAL *8(a-h,o-z)
	    real *8 p(2,26)
        integer n
        integer i,jj
        real *8 rtek,z
        
	    do i=2,n
	        if (p(1,i).ge.z) goto 1
	    end do
	    i=n
 1  	continue
	    rtek=p(2,i-1)+(p(2,i)-p(2,i-1))/(p(1,i)-p(1,i-1))*(z-p(1,i-1))
	    return
	end subroutine interp
	
!    function defFactor(Re,u,nu)
!    implicit none
!        real *8 defFactor
!        real *8 Re,u,nu,xlam
!        if (Re.le.1500d0) then
!            defFactor=2.0d0
!        else
!            call def_xlam(xlam,u,Re,nu)
!            defFactor=1d0/xlam
!        end if
!        return
!    end function defFactor
        
    subroutine def_xlam(xlam,u,r,mu)
        real *8 xlam,u,r,mu
        real *8 x1,x,x2,f
    
        x1=1d-5
        x2=1.0
        n=0
    1	x=(x1+x2)/2.0d0
        f=1.0d0/dsqrt(x)-2.0d0*dlog10(u*2.0d0*r/mu*dsqrt(x))+0.8d0
        if (f.gt.0) x1=x
        if (f.lt.0) x2=x
        n=n+1
        if (n.lt.100) goto 1
        xlam=x
    return
    end  subroutine def_xlam
    
    subroutine add_point(n,l,m,k,nodeInterface)
    implicit none
        integer,intent(in) :: l,m,nodeInterface
        integer,intent(inout) :: n(nodeInterface)
        integer,intent(inout) :: k
        
        integer :: j

!        if (k.gt.0) then
!            do j=l+1,l+k
!                if (n(j).eq.m) return
!            end do
!        end if
!        k=k+1
!        n(l+k)=m
        
        
        do j = 1+l,nodeInterface+l
            if (n(j).eq.m) return
            if(n(j).eq.0)then
                k = k + 1
                n(j) = m
                return
            end if
        end do
        
        !	write(*,*) 'n(l+k),k,l=',n(l+k),k,l
        return
    end subroutine add_point
    
        !	call renumb(nodeBoundaryMeshLabel,pointProcAddressing,nodeInterfacet,m,nPoints)
    subroutine renumb(n,m,l,k,j1)
        !nodeBoundaryMeshLabel,pointProcAddressing,nodeInterfacet,m,nPoints
    implicit none
        integer,intent(in) :: m(1),l,k,j1
        integer,intent(inout) :: n(j1)
        
        integer :: i,j,ni
        
        do i=l+1,l+k-1
            do j=i+1,l+k
                if (n(i).gt.n(j)) then
                    ni=n(i)
                    n(i)=n(j)
                    n(j)=ni
                end if
            end do
        end do
!        do i=l+1,l+k
!            do j=1,j1
!                if (m(j).eq.n(i)) n(i)=j
!            end do
!        end do
        
        do i=l+1,l+k
            do j=1,j1
                if (m(j).eq.n(i))then
                    n(i)=j
                    exit
                end if
            end do
        end do
        
    return
    end subroutine renumb

    subroutine intervt(p,n,z,v,t,pres)
    implicit none
        !IMPLICIT REAL *8(a-h,o-z)
        real(8),intent(in) :: p(4,16)
        integer,intent(in) :: n
        real(8),intent(in) :: z
        real(8),intent(out) :: v,t,pres
        
        integer i
        
        do i=2,n
            if (p(1,i).ge.z) goto 1
        end do
        i=n
    1	continue
        v=p(2,i-1)+(p(2,i)-p(2,i-1))/(p(1,i)-p(1,i-1))*(z-p(1,i-1))
        t=p(3,i-1)+(p(3,i)-p(3,i-1))/(p(1,i)-p(1,i-1))*(z-p(1,i-1))
		pres=p(4,i-1)+(p(4,i)-p(4,i-1))/(p(1,i)-p(1,i-1))*(z-p(1,i-1))

    return
    end subroutine intervt

    subroutine intervt_old(p,n,z,v,t)
    implicit none
        !IMPLICIT REAL *8(a-h,o-z)
        real(8),intent(in) :: p(3,16)
        integer,intent(in) :: n
        real(8),intent(in) :: z
        real(8),intent(out) :: v,t
        
        integer i
        
        do i=2,n
            if (p(1,i).ge.z) goto 1
        end do
        i=n
    1	continue
        v=p(2,i-1)+(p(2,i)-p(2,i-1))/(p(1,i)-p(1,i-1))*(z-p(1,i-1))
        t=p(3,i-1)+(p(3,i)-p(3,i-1))/(p(1,i)-p(1,i-1))*(z-p(1,i-1))
    return
    end subroutine intervt_old

    subroutine GetAlfa(x, alfa)
    implicit none
        real(8),intent(in) :: x
        real(8),intent(inout) :: alfa
        real(8) :: epsilon, halfPi
        
        epsilon = 1.0d-7
        
        if(x.lt.-epsilon)then
            alfa = 1.0d0
            return
        end if
        if(x.gt.epsilon)then
            alfa = 0.0d0
            return
        end if

!        alfa = 0.5 * (1.0d0 - x/epsilon)
!        if(x.lt.0) alfa = 1.0d0
        halfPi = 1.5707963267948966192313216916398
        alfa = 0.5*(-sin(halfPi*x/epsilon) + 1.0d0)
        return
    end subroutine


subroutine CalcBoundaryFaceCubes()
use Data_cfoam
use DataTime
use Data_Mesh
implicit none

    integer :: iCell,iCellFace, iFace, i, OpenFile, n
    real(8) :: Distance, r(3), Rxy
    logical :: flag

    pnew = -999.0d0
    Unew = -999.0d0
    Tsnew = -999.0d0

    ! Nozzle 1
    do iFace = nInternalFaces+1,nFaces
!        r(1:3) = meshSf(1:3,iFace)
!        if(abs(abs(r(1))-1d0).lt.1d-5)then
            r(1:3) = meshFaceCenters(1:3,iFace)

        if(abs(r(1)-0.6d0).lt.1d-4)then

!        if(abs(r(1)+0.6d0).lt.5d-4)then

            Rxy = sqrt(r(2)**2+r(3)**2)
!            if(Rxy-0.58.lt.0)then
!            if(Rxy-0.604d0.lt.0)then
!            if(Rxy-1.4d0.lt.0)then
            if(Rxy-0.6405d0.lt.0)then
                iCell = meshFaceOwner(iFace)
                pnew(iCell) = iFace
                Unew(1:3,iCell) = r(1:3)
                Tsnew(iFace) = 1d0
            end if
        end if
    end do
    ttime = ttime + 1d0
    if(OutputFromOneProc)then
        call OutputFromAllProcs()
    else
        call WriteRes()
    end if
    ! ******************************************* !

    pnew = -999.0d0
    Unew = -999.0d0

    ! Bypass 2
    do iFace = nInternalFaces+1,nFaces
        if(Tsnew(iFace).gt.0)then
!            write(*,*) '@warning: Bypass'
            cycle
        end if

        r(1:3) = meshFaceCenters(1:3,iFace)

!        if(abs(r(1)+1.5d0).lt.1d-5)then
        if(abs(r(1)+0.3d0).lt.5d-4)then
        

!        r(1:3) = meshSf(1:3,iFace)
!        if(abs(abs(r(1))-1d0).lt.1d-5)then
            r(1:3) = meshFaceCenters(1:3,iFace)

            Rxy = sqrt(r(2)**2+r(3)**2)
!            if(Rxy.gt.0.693.and.Rxy.lt.1.29)then
!            if(Rxy.gt.0.79.and.Rxy.lt.1.076)then
            if(Rxy-1.2d0.lt.0)then
                iCell = meshFaceOwner(iFace)
                pnew(iCell) = iFace
                Unew(1:3,iCell) = r(1:3)
                Tsnew(iFace) = 2d0
            end if
        end if
    end do
    ttime = ttime + 1d0
    if(OutputFromOneProc)then
        call OutputFromAllProcs()
    else
        call WriteRes()
    end if
    ! ******************************************* !

    pnew = -999.0d0
    Unew = -999.0d0

    ! inlet 3
    do iFace = nInternalFaces+1,nFaces
        if(Tsnew(iFace).gt.0)then
!            write(*,*) '@warning: inlet'
            cycle
        end if
        r(1:3) = meshFaceCenters(1:3,iFace)

!        if(r(1).lt.-3.8d0)then
        if(r(1).lt.-2.8d0)then
            iCell = meshFaceOwner(iFace)
            pnew(iCell) = iFace
            Unew(1:3,iCell) = r(1:3)
            Tsnew(iFace) = 3d0
        end if
    end do
    ttime = ttime + 1d0
    if(OutputFromOneProc)then
        call OutputFromAllProcs()
    else
        call WriteRes()
    end if

    ! ******************************************* !

    pnew = -999.0d0
    Unew = -999.0d0

    ! no-slip BC 4
    do iFace = nInternalFaces+1,nFaces
        if(Tsnew(iFace).gt.0) cycle
        r(1:3) = meshFaceCenters(1:3,iFace)
        Rxy = sqrt(r(2)**2+r(3)**2)

!        if(r(1).lt.10d0)then
!            if(dabs(r(3)).lt.14.5d0)then
!                if(dabs(r(2)).lt.2d0)then
        if(r(1).lt.5d0)then
            if(dabs(r(3)).lt.3d0)then
                if(dabs(r(2)).lt.3d0)then
                    iCell = meshFaceOwner(iFace)
                    pnew(iCell) = iFace
                    Unew(1:3,iCell) = r(1:3)
                    Tsnew(iFace) = 4d0
                end if
            end if
        end if
    end do
    ttime = ttime + 1d0
    if(OutputFromOneProc)then
        call OutputFromAllProcs()
    else
        call WriteRes()
    end if
    ! ******************************************* !

    pnew = -999.0d0
    Unew = -999.0d0

    ! outlet BC 5
    do iFace = nInternalFaces+1,nFaces
        if(Tsnew(iFace).lt.0d0)then
            iCell = meshFaceOwner(iFace)
            pnew(iCell) = iFace
            Unew(1:3,iCell) = r(1:3)
            Tsnew(iFace) = 5d0
        end if
    end do
    ttime = ttime + 1d0
    if(OutputFromOneProc)then
        call OutputFromAllProcs()
    else
        call WriteRes()
    end if
    ! ******************************************* !

    pnew = -999.0d0
    Unew = -999.0d0

    ! all boundary faces 6
    do iFace = nInternalFaces+1,nFaces
        iCell = meshFaceOwner(iFace)
        pnew(iCell) = iFace
        Unew(1:3,iCell) = r(1:3)
    end do
    ttime = ttime + 1d0
    if(OutputFromOneProc)then
        call OutputFromAllProcs()
    else
        call WriteRes()
    end if
    ! ******************************************* !


! write new face and owner according to new order of boundary faces

    ! owner_new
        open(2,file='owner')
        write(2,*) ' ' !/*--------------------------------*- C++ -*----------------------------------*\'
        write(2,*) '| =========                 |                                                 |'
        write(2,*) '| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |'
        write(2,*) '|  \\    /   O peration     | Version:  1.7.1                                 |'
        write(2,*) '|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |'
        write(2,*) '|    \\/     M anipulation  |                                                 |'
        write(2,*) '\*---------------------------------------------------------------------------*/'
        write(2,*) 'FoamFile'
        write(2,*) '{'
        write(2,*) '    version     2.0;'
        write(2,*) '    format      ascii;'
        write(2,*) '    class       labelList;'
        write(2,*) '    note        "nPoints:',nPoints, ' nCells:', nCells, &
                    ' nFaces:',nFaces, ' nInternalFaces:',nInternalFaces, '";'
        write(2,*) '    location    "constant/polyMesh";'
        write(2,*) '    object      owner;'
        write(2,*) '}'
        write(2,*) '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
        write(2,*) ' '
        write(2,*) ' '
        write(2,*) nFaces
        write(2,*) '('

        do iFace = 1,nInternalFaces
!            write(2,*) meshFaceOwner(iFace2iFaceNew(iFace))-1
            write(2,*) meshFaceOwner(iFace)-1
        end do
        
        do i = 1,5
            do iFace = nInternalFaces+1,nFaces
                if(abs(Tsnew(iFace)-real(i)).lt.1d-1)then
                    write(2,*) meshFaceOwner(iFace)-1
                end if
            end do
        end do

        write(2,*) ')'
        write(2,*) ' '
        write(2,*) ' '
        write(2,*) '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'

        close(2)

! ***************************************************************** !

! faces_new
        open(2,file='faces')
        write(2,*) ' ' !/*--------------------------------*- C++ -*----------------------------------*\'
        write(2,*) '| =========                 |                                                 |'
        write(2,*) '| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |'
        write(2,*) '|  \\    /   O peration     | Version:  1.7.1                                 |'
        write(2,*) '|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |'
        write(2,*) '|    \\/     M anipulation  |                                                 |'
        write(2,*) '\*---------------------------------------------------------------------------*/'
        write(2,*) 'FoamFile'
        write(2,*) '{'
        write(2,*) '    version     2.0;'
        write(2,*) '    format      ascii;'
        write(2,*) '    class       faceList;'
        write(2,*) '    location    "constant/polyMesh";'
        write(2,*) '    object      faces;'
        write(2,*) '}'
        write(2,*) '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
        write(2,*) ' '
        write(2,*) ' '
        write(2,*) nFaces
        write(2,*) '('
        
        do iFace = 1,nInternalFaces
!            write(2,*) '4(',meshFaces(1:4,iFace2iFaceNew(iFace))-1,')'
            write(2,*) '4(',(meshFaces(1:4,iFace))-1,')'
        end do

        do i = 1,5
            do iFace = nInternalFaces+1,nFaces
                if(abs(Tsnew(iFace)-real(i)).lt.1d-1)then
                    write(2,*) '4(',(meshFaces(1:4,iFace))-1,')'
                end if
            end do
        end do

        write(2,*) ')'
        write(2,*) ' '
        write(2,*) ' '
        write(2,*) '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'

        close(2)

! ***************************************************************** !

    ! output face numbers for every boundary
    open(11,file='BoundaryFaceNumbers.dat',IOSTAT=OpenFile)
    if (OpenFile.ne.0) then
	    write(*,*) '@@ Fatal Error (CalcBoundaryFaceCubes): File open error '//processor(1:ie)//'ResCells'//String7Number(1:7)
	    STOP '@@ Fatal Error!'
    end if

    write(11,*) 'nBoundaryFaces=', nFaces - nInternalFaces
    write(*,*) 'nBoundaryFaces=', nFaces - nInternalFaces
    do i = 1,5
        write(11,*) 'Boundary ', i

        flag = .false.
        n = 0
        do iFace = nInternalFaces+1,nFaces
            if(abs(Tsnew(iFace)-i).lt.1d-10)then
                n = n + 1
                if(flag) cycle
                flag = .true.
                write(11,*) iFace-1
            else
                if(flag)then
                    flag = .false.
                    write(11,*) iFace-1
        !            write(11,*) '------------'
                end if
                
                if(Tsnew(iFace).lt.0)then
                    write(*,*) '@warning: Tsnew(iFace)<0, iFace=', iFace
                end if
            end if
        end do
        if(flag)then
            write(11,*) iFace-1
        end if
        
        write(11,*) 'nFace=',n
        write(*,*) 'nFace=',n
        
        write(11,*) '***************************************************************'
    end do
    
    close(11)
    
!!    pnew = -999.0d0
!!    Unew = -999.0d0

!!    ! all boundary faces
!!    do iFace = nInternalFaces+1,nFaces
!!        iCell = meshFaceOwner(iFace)
!!        pnew(iCell) = iFace
!!        Unew(1:3,iCell) = r(1:3)
!!        Tnew(iCell) = 6
!!    end do
!!    ttime = ttime + 1d0
!!    if(OutputFromOneProc)then
!!        call OutputFromAllProcs()
!!    else
!!        call WriteRes()
!!    end if

    return

! ****************************************************************
    ! outlet
        do iCell = 1,nCells
            do iCellFace = 1,6
                iFace = meshCells(iCellFace,iCell)
                if(iFace.gt.nInternalFaces)then
                    if(meshFaceCenters(1,iFace).gt.480d0)then
                        pnew(iCell) = iFace
                        Unew(1:3,iCell) = meshFaceCenters(1:3,iFace)
                    end if
                end if
            end do
        end do
        ttime = ttime + 1d0
        if(OutputFromOneProc)then
            call OutputFromAllProcs()
        else
            call WriteRes()
        end if
    ! ******************************************* !
    ! walls
        pnew = -999.0d0
        Unew = -999.0d0
        do iCell = 1,nCells
            do iCellFace = 1,6
                iFace = meshCells(iCellFace,iCell)
                if(iFace.gt.nInternalFaces)then
                    if(meshFaceCenters(2,iFace).lt.-0.9d0)then
                        pnew(iCell) = iFace
                        Unew(1:3,iCell) = meshFaceCenters(1:3,iFace)
                    end if
                end if
            end do
        end do
        ttime = ttime + 1d0
        if(OutputFromOneProc)then
            call OutputFromAllProcs()
        else
            call WriteRes()
        end if

        pnew = -999.0d0
        Unew = -999.0d0
        do iCell = 1,nCells
            do iCellFace = 1,6
                iFace = meshCells(iCellFace,iCell)
                if(iFace.gt.nInternalFaces)then
                    if(meshFaceCenters(2,iFace).gt.-0.1d0)then
                        pnew(iCell) = iFace
                        Unew(1:3,iCell) = meshFaceCenters(1:3,iFace)
                    end if
                end if
            end do
        end do
        ttime = ttime + 1d0
        if(OutputFromOneProc)then
            call OutputFromAllProcs()
        else
            call WriteRes()
        end if

    ! ******************************************* !
    ! periodic
        pnew = -999.0d0
        Unew = -999.0d0
        do iCell = 1,nCells
            do iCellFace = 1,6
                iFace = meshCells(iCellFace,iCell)
                if(iFace.gt.nInternalFaces)then
                    if(meshFaceCenters(1,iFace).lt.480d0.and. &
                        meshFaceCenters(1,iFace).gt.1d0.and. &
                        abs(meshFaceCenters(3,iFace)).lt.0.1d0.and. &
                        meshFaceCenters(2,iFace).gt.-0.9d0.and. &
                        meshFaceCenters(2,iFace).lt.-0.1d0)then
                            pnew(iCell) = iFace
                            Unew(1:3,iCell) = meshFaceCenters(1:3,iFace)
                    end if
                end if
            end do
        end do
        ttime = ttime + 1d0
        if(OutputFromOneProc)then
            call OutputFromAllProcs()
        else
            call WriteRes()
        end if
    ! ******************************************* !
    ! airfoil
        pnew = -999.0d0
        Unew = -999.0d0
        do iCell = 1,nCells
            do iCellFace = 1,6
                iFace = meshCells(iCellFace,iCell)
                if(iFace.gt.nInternalFaces)then
                    if(meshFaceCenters(1,iFace).lt.1d0.and. &
                        meshFaceCenters(1,iFace).gt.-1d0.and. &
                        abs(meshFaceCenters(3,iFace)).lt.10d0.and. &
                        meshFaceCenters(2,iFace).gt.-0.9d0.and. &
                        meshFaceCenters(2,iFace).lt.-0.1d0)then
                            pnew(iCell) = iFace
                            Unew(1:3,iCell) = meshFaceCenters(1:3,iFace)
                    end if
                end if
            end do
        end do
        ttime = ttime + 1d0
        if(OutputFromOneProc)then
            call OutputFromAllProcs()
        else
            call WriteRes()
        end if
    ! ******************************************* !
    ! inlet
        pnew = -999.0d0
        Unew = -999.0d0
        do iCell = 1,nCells
            do iCellFace = 1,6
                iFace = meshCells(iCellFace,iCell)
                if(iFace.gt.nInternalFaces)then
                    Distance = sqrt(sum(meshFaceCenters(1:3,iFace)*meshFaceCenters(1:3,iFace)))
                    if((meshFaceCenters(1,iFace).lt.-1d0.or.abs(meshFaceCenters(3,iFace)).gt.100d0).and. &
                        meshFaceCenters(2,iFace).gt.-0.9d0.and. &
                        meshFaceCenters(2,iFace).lt.-0.1d0.and. &
                        meshFaceCenters(1,iFace).lt.480d0)then
                            pnew(iCell) = iFace
                            Unew(1:3,iCell) = meshFaceCenters(1:3,iFace)
                    end if
                end if
            end do
        end do
        ttime = ttime + 1d0
        if(OutputFromOneProc)then
            call OutputFromAllProcs()
        else
            call WriteRes()
        end if
    ! ******************************************* !
    
!        do iCellFace = 1,6
!            do iCell = 1,nCells
!                iFace = meshCells(iCellFace,iCell)
!                if(iFace.gt.nInternalFaces)then
!                    if(iCellFace.eq.1)then
!                        p(iCell) = iFace
!                        U(1:3,iCell) = meshFaceCenters(1:3,iFace)
!                    end if
!                else
!                    p(iCell) = -999.0d0
!                    U(1:3,iCell) = -999.0d0
!                end if
!            end do
!            ttime = ttime + 1d0
!            call OutputFromAllProcs()

!            write(32) ttime,0,nCells !,(U(1:3,i),p(i),T(i),i=1,nCells) !,E(i),rho(i)
!            write(32) (U(1:3,i),p(i),T(i),i=1,nCells) !,E(i),rho(i)
!        end do
        Stop 'CalcBoundaryFaceCubes'
end subroutine CalcBoundaryFaceCubes

subroutine OutputFaceCubes()
use Data_cfoam
use DataTime
use Data_Mesh
implicit none

    integer :: j, iFace, iCell

! processor number output
    p = rank
    U = rank
    
    ttime = ttime + 1d0
    
    pnew = p
    Unew = U
    Tnew = 0d0       
    
    if(OutputFromOneProc)then
        call OutputFromAllProcs()
    else
        call WriteRes()
    end if

    do j=1,nBoundaryFields
        p = -999.0d0
        U = -999.0d0
        do iFace=1+boundaryMeshStart(j),boundaryMeshStart(j)+boundaryMeshSize(j)
            iCell = meshFaceOwner(iFace)
            p(iCell) = iFace
            U(1:3,iCell) = meshFaceCenters(1:3,iFace)
        end do
        ttime = ttime + 1d0
        
        pnew = p
        Unew = U
        Tnew = 0d0       
        
        if(OutputFromOneProc)then
            call OutputFromAllProcs()
        else
            call WriteRes()
        end if
    end do


    Stop 'OutputFaceCubes'
        
end subroutine OutputFaceCubes

subroutine StopSimulation()
use TimeMeasure
implicit none

    call outputHWTimer()
    Stop 'SignalC works good!'

end subroutine StopSimulation
